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Abstract 

The identihcation of relevant cohective coordinates is crucial for the interpretation of coherent 
nonlinear spectroscopies of complex molecules and liquids. Using an h expansion of Liouville space 
generating functions, we show how to factorize multitime nonlinear response functions into prod- 
ucts of lower-order correlation functions of collective coordinates, and derive closed expressions for 
linear, second and third order response functions. In addition to providing systematic quantum 
corrections, h offers a convenient bookkeeping device even for the purely classical response, since 
including quantum fluctuations allows to circumvent the expensive computation of stability matri- 
ces which is a major bottleneck in Molecular Dynamics (MD) simulations. The existing classical 
simulation strategies, including Mode-Coupling in k space and in real-space, Langevin equations, 
and Instantaneous Normal Modes are compared from a unified viewpoint. 
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I. INTRODUCTION 



In 1993 Tanimura and Mukamel had proposed the fifth order Raman response^ as a 
multidimensional spectroscopic technique especially suitable for investigating the structure 
and dynamics of molecular liquids by revealing detailed information unavailable from linear 
spectroscopies. That article had triggered an intense experimental investigations mainly on 
liquid CS2- 2,3,4,5,6,7,8,9,10,11,12,13,14,15,16 Eariigj; experimental investigations were haunted by 

competing, sequential, low order, (cascading) processes. Separating the direct and sequential 
contributions had drawn considerable attention. ^'^^'^^'^^ 

Several theoretical methods have been employed to predict the fifth order response from 
molecular dynamics simulations of liquids. Two methods obtain the response directly with- 
out further approximations, other than that the response is classical. The first, based on 
calculating time correlation functions, relies on propagation of the full stability matrix (Eq. 
(25)).^^'^°'^^ Since the stability matrix depends on the number of phase space coordinates 
squared, it is very time consuming and was only implemented for very small systems^^'^^'^^'^^. 
The other. Finite Field, method is based on propagating only one column of the stabihty 
matrix, giving rise to a particular response function of interest, significantly reducing the 
computational effort. This method is a direct simulation of the experiment, where forces 
originating from interactions between the electric fields and the molecules are incorporated 
in the simulation on the fly. 26,27,28 Qj^^ drawback of computing the actual non equilibrium 
response rather than response functions is that the entire the full simulation needs to be 
repeated for each choice of time intervals and pulse configurations. Both of these methods 
are therefore computationally very demanding. Developing alternative numerically more af- 
fordable approaches, which could provide physical insight will therefore be highly desirable. 

Some discrepancies currently exist between various real-space simulations performed un- 
der slightly different simulation conditions on liquid C5'2.^^'^^ These differences are most 
pronounced along the second time axis, (T32), where both a ridge^^'^° and nodes^^'^^ have 
been reported. A fundamental understanding of the underlying physical processes should 
help resolve the questions about the origin of the nodes and the ridges. 

The first approximate scheme employed to analyze the fifth order Raman response was 
based on the Instantaneous Normal Modes (INM) 20,22,31,32,33,34,35,36,37^ This method uses 
snapshots of the liquid "normal modes" assuming that they are harmonic and do not change 
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over the timescale of the experiment. However, in general the normal modes do change on 
the timescale of their own periods^^ and recent studies have shown that the INM gives rather 
poor results for the fifth order response of liquid Xenon .^^'^^ 

A more attractive procedure is to identify some relevant collective coordinates and adopt 
a reduced description for the response. Unlike microscopic INM where a harmonic model 
for molecular liquids may not be justified, collective coordinates can have Gaussian statis- 
tics by virtue of the central limit theorem. A simple and tractable physical picture for the 
origin of the response is then, in principle, possible. The multimode Brownian oscillator 
model has been successfully employed in the analysis of solvation dynamics in electronic 
spectroscopy,^^ where the response may be expressed using a few (overdamped or 
underdamped) collective coordinates. This model has been used to simulate the fifth or- 
der response of liquid water but identifying the microscopic origin of these modes still 
remains an open challenge^^. Nonlinear hydrodynamics and mode coupling theories suc- 
cessfully use collective variables in momentum (k) space to describe slow, long- wavelength, 
variables and their fluctuations. Mode-Coupling (MC) theory^^'^^'^'^'^^'^^'^'^'^-'^ has been ap- 
plied to relate the fifth order Raman response^^'^^'^^'^^'^^ to fluctuations of density modes is k 
space. Another related approach is based on the Generalized Langevin Equations (GLE)^^'^^. 
Classical mode-coupling theory contains some ambiguities regarding the proper factoriza- 
tion of high order response functions, and Schofield^^ and Keyes^^ and collaborators have 
discussed possible simulation strategies based on Langevin equations. 

In this paper we apply a unifying picture of quantum field and mode-coupling Green 
function theories developed recently^^'^*^'^^ to derive expressions for the first-, second- and 
third order response functions. The technique provides an unambiguous and unique factor- 
ization scheme of multitime correlation functions and allows the perturbative incorporation 
of anharmonicities as well as quantum corrections through an h expansion. Applications are 
made to the fifth order Raman response and compared with other approximate methods. 
In Section II we present the superoperator formalism. In Section III we describe how third 
and fifth order Raman response can be obtained from the general first- and second order 
response functions. We discuss the connections, similarities and differences of the present 
formulation with other approaches. Conclusions are given in Section IV. 



3 



II. LIOUVILLE SPACE FORMULATION OF RESPONSE FUNCTIONS 



The present approach is based on superoperators (^4+ and A^) corresponding to an 
ordinary Hilbert space operator A, defined by their action on some Hilbert space operator 



A+n = -{AQ + QA), 

A_n = An- nA. (i) 

Using this notation, A compact expression for the n dimensional (nD) quantum mechani- 
cal response functions can then be written in terms of the dipole superoperators and 

^ ^ 60,61,62 

^^"""^ = (^)" WK+l)/^-(rn) ■ ■ -/^-(n)) . (2) 

The superscript of i?^"^) indicates that it depends on the n time intervals between successive 
times Ti ■ ■ ■ r„_|_i and thus constitutes an n dimensional technique. The average (A) = 
Tr[Apeq] denotes the trace with respect to the equilibrium density matrix of the system p^g. 

The Hamiltonian will be partitioned into a Harmonic, quadratic, part {Hq) and an an- 
harmonic (V) part 

H^Ho + V, (3) 

and the response function of a weakly anharmonic system will be expanded perturbatively 
in V.^^ 

°° / 7 \ /■*"+! 



h 

m=0 ^ / -"-^ 

X (rMT„H-i)/i_(T„)---/i4Ti)y_(T;)--.T/_(T0)^. (4) 

Here Vi, and jlv are the superoperators associated with the anharmonic part of the potential 
and with the interaction dipole, respectively in the interaction picture with respect to Hq, 
i.e., 

A^{t) = exp (^LoTj Ai^exp (-^LotJ . (5) 
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Here 

LoO = [Ho,n] = {Ho)-n, (6) 

is the Liouville operator corresponding to Hq . The average {A)o = Tr[Apo] is defined as 
the trace with respect to the equihbrium density operator for the harmonic system po- T is 
the time ordering operator in Liouville space which arranges all superoperators so that their 
time arguments decrease from left to right^°. 

We shall represent Hq in terms a few primary (collective) coordinates Qj'^^'^^'^^ described 
by the Hamiltonian 

where Pj{Qj) is the momentum (coordinate) operator of the j'th primary mode, flj and Mj 
are its frequency and reduced mass respectively. The anharmonic potential V, is 

nQ) = j:j^K''-lQn---Q3.- (8) 

N=3 

We assume that the dipole operator /i only depends on the primary coordinates and expand 
it as 

oo ^ 

^'=T.J^^n-LQn---Qj.- (9) 

N=l 

The primary modes further interact with a large number of low- frequency harmonic (bath) 
coordinates which induce relaxation and dephasing. These bath degrees of freedom q and 
their linear coupling to the primary modes are described by the Hamiltonian Hb- 
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and PjaiQja) are momentum (coordinate) operators of bath oscillators. Cja are the coupling 
constants between the primary and bath coordinates. The total harmonic Hamiltonian is 

, 40 

given by 

//o = ^^m(Q) + i^B(Q,q). (11) 

Applying the algebraic rules for the superoperators given in Appendix A^° the Taylor 
expansions of //+(t), //-(r) and T^-(t) can be expressed in terms of the elementary superop- 
erators Qjj^ and Qj-- Using these expansions, we can convert the time-ordered product of 



superoperators into a time ordered product of primary superoperators. This transforms the 
computation of response functions (Eq. (4)) to the evaluation of products of the form 

Wij^l^rnTm} = {QjivA^l) ■ ■ ■ Qjjvi^jv ('^Jv))o, (12) 

where i^i^...^ i^n — i, and jm runs over the collective coordinates. Note that the number N of 
operators in the product needed to compute /?("^) is generally greater than n,N > n. The 
reasons are (i) jj-i, may be nonlinear in the elementary operators, (ii) The expansion in V- 
adds more operators to the product. Generally some of the tj in Eq. (12) will be the same 
since i?^"^^) only depends on n + m + 1 times, which is smaller than or equal to N. 
We next introduce the superoperator generating functional^^ 



S{{J{t)})=(Tcxp 



Yl / Jju{r)Qju{r)dr 



(13) 



Since the Hamiltonian Hq is quadratic, the generating functional may be computed ex- 
actly using the second order cumulant expansion. This gives 

{/•CO l>T2 
J2 drj dr, (14) 
j,k -'-'^ "'-'^ 

[-ihJ,+ {T-,)Jk-{T,) G+-(r2i) + J,+(T2)4+(ri)G++(T2i) 

We have introduced the notation = — Tj and the two basic Liouville space Green 
functions. 

Gti.r,,) ^ '-{TQX{T,)Ql{n))^, (15) 

GS+(r,i) ^ {TQX{r,)Qi{n))^. (16) 

Using the Eq. (1), Eqs.(15) and (16) can be recast as combinations of ordinary (Hilbert 
space) correlation functions 

G^-{t,,) = e{r,^)'-{{Q\T,)Q^{r,))^ - {Q^rM {r,)) ^) , (17) 

Gttir,,) = \{{Q\r,)Q={n)), + {QKri)Q\r,))^). (18) 

and may also be expressed in terms of the spectral densities Cij{uj) 

/OO J 
J^CMM^rl (19) 



OO 



Gt^{r) = h I ^Qjiu) cos(u;r) coth ( ) . (20) 



Eq. (19) is the definition of the spectral density, where 9{t) is the Heaviside step function 
(equal to zero for r < and equal to one for r > 0) . In the classical, high temperature, 
hmit coth.{huj/2kBT) ^ 2kBT/hw, and both and become independent of h. In 
that case the two are related by the classical fluctuation-dissipation relation 

Gt(r) = -^^(r)±GtHr). (21) 

Semiclassical approximations to the response may be developed by expanding G++ in powers 
of h. 

Time-ordered correlation functions of superoperators may be obtained from the generat- 
ing functional by functional derivatives^^ 

In order to compute the response function (which gives the response to very short pulses) 
to a given order in the field, the generating functional can be simplified since only a limited 
number of times will contribute, and the primary operators connected with the last time 
will have to be -I- operators since {A_{ti)B^{t2)) vanishes for ti ^ T2. ^° The generating 
functional for two time quantities {Ai,{ti)B^i{t2)) thus reads 

S^'\{J{r)}) = exp { ^ -iH-+(T2) Jfe-(Tl)G}r(T2l) + J, + (t2) Jfe+(Ti)G++(T2l) 

+ J,+(t2) Jfe+(T2)G++(0) + J,-+(ti) Jfe+(Ti)G++(0)}. (23) 

Here wc assumed Gaussian statistics of Qj so that the exact generating functional is given 
by the second order cumulant expansion. The generating functional provides a compact 
form for Wick's theorem. Generating functionals for multitime quantities may be written in 
a similar way. 

Using the general expression for the n'th order response function given in Appendix A, 
we can expand i?^"^) to any desired order in the primary operators. Since the J's always 
come in pairs in the generating functional, the derivatives will vanish for all terms with an 
odd number of elementary operators, once the J = limit is taken. 

Note the delicate interplay of the h factors in Eq. (A6). Keeping h alive even in the 
classical limit, is what allows us to avoid the computation of stability matrices, h retains 
information about quantum fluctuations which are differences between "left" and "right" 
trajectories^^'^^'^^: The stabihty matrices are the corresponding derivatives as h tends to 



zero. In the classical evaluation of Eq. (A6) we keep terms in the generating function to order 
^^n+m'^ ^ ^jjgj^ cancels out by the prefactor in Eq. (4) and the result is independent of h, as it 
should be.^^ Higher order terms in h provide quantum corrections to the response. We further 
note that the response function, which is a particular combination of correlation functions, 
has a well-defined classical limit. Individual correlation functions in Hilbert space generally 
do not have a clear physical meaning and consequently their h ^ limit is ill defined. H 
therefore cancels only once these combinations are evaluated to yield an observable. 

As an example, the ID response function expanded to fourth order in the elementary 
operators is given by 

ij 

+ E4'^/^i?^.V(^2i)G++(r,0 + E/^SW'^^^(^2i)^i^(0)- (24) 

ijkl ijkl 

These two lowest-orders in the expansion are independent of the anharmonicity and the ID 
response can be expected to be dominated by the harmonic part of the potential. Closed 
expressions for the 2D and 3D response functions expanded to sixth order in the primary 
coordinates are given in Appendix B. 

III. COMPARISON OF SIMULATION STRATEGIES FOR 2D FIFTH ORDER 

RAMAN RESPONSE 

In off-resonant Raman spectroscopy, the field-matter interaction comes through the dipole 
moment induced by an electric field instead of the permanent dipole moment. The Raman 
response can therefore be obtained by simply substituting the dipole operators in the ex- 
pressions in Eq. (24) and Appendix B with the operator for the induced dipole, which in 
turn is given by the product of the polarizability and the inducing field, /i is therefore simply 
replaced hy a- E. This substitution leads to a higher order dependence on the electric field: 
The nD response is n'th order in the field for dipole (e.g. infrared) response but (2n + l)'th 
order for Raman. The Raman response that is third order in the electric fields is therefore 
described by the ID response function. The fifth order Raman response is given by the 2D 
response function and so forth. 

The various approaches used for the simulation of fifth-order Raman signals differ not 
only by the simulation technique, but also by the model used, which complicates their direct 
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comparison. The expressions derived in this paper allow a critical comparison of the various 
simulations. This will be done next. 

Classical MD simulations in real-space: 

The real-space simulations by Saito et al?*^ Ma et al?"^'"^^ and Jansen et a/.^^'^'' use 
the nuclear coordinates as basis and the response is calculated without invoking Wick's 
theorem. The calculated classical response functions therefore formally include all orders 
of the expansion derived here. The price is the need to compute stability matrices. The 
nuclear coordinates might not be the best choice of basis for the expansion derived in the 
previous section and analyzing the response in this basis might be inconvenient. The real- 
space simulations give the full response including all anharmonicities and nonlinearitics of 
the polarizability. Analyzing the response calculated with these methods is a daunting task, 
since all terms which depend on different modes and nonlinearitics are added. This makes 
it difficult to establish the connection between the spectral features and the underlying 
dynamics. 

The stability matrix (M(t2,ti)) is a A?" by matrix, where A" is the number of phase 
space coordinates. ^^'^'^'^^'^^ Each matrix element is given by the derivative of a phase space 
coordinate at time ri with respect to a phase space coordinate Xj at another time, T2 



Schemes for obtaining the stability matrix using the Hessian matrix have been described in 



Equilibrium simulations 19.20,22,23,24,25 based on propagating the full stabihty matrix^^ 
in order to evaluate the Poisson bracket arising in the time correlation function expression for 
the fifth-order response. The major bottleneck in the equilibrium method is the propagation 
of the N xN stability matrix. ^^'^'^ In contrast, in the non equilibrium. Finite Field, approach 
26,27,64^ propagation of the full stability matrix is avoided and only a few vectors of dimension 
A^ corresponding to the response function are propagated. In this approach the evaluation 
of the first order derivatives of the polarizabihty, needed in order to calculate the forces 
exerted by the electric fields, is the most time consuming part of the simulation. 

Instantaneous Normal Modes: 

For a collection of oscillators with frequencies coi we have^ 




(25) 



the literature. i9'20'23>25,64,65,66 



1 



MiUJi 
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1 



cos(a;jr) coth 



( 



2hT 



hhJi 




(26) 



INM simulations combine these expressions with MD simulations of real hquids. The system 
is assumed to evolve independently in the different normal modes giving rise to the Kronecker 
deltas in Eq. (26). INM simulations have most often been applied in the classical, high 
temperature limit, where only the first term in the expansion of coth is retained (coth(a;) ~ 
1/x). Eqs. (B1)-(B6) reduce to the INM expressions if the INM coordinates are used as the 
primary operators and the anharmonicities are neglected. In this paper we use a different 
bookkeeping: terms have been kept to a certain order in the primary operators, whereas 
the INM expressions traditionally have been truncated by neglecting terms containing third- 
and higher order nonlinearities of the polarizabilty 20,22,36,37^ retaining only Eqs. (B2) and 
(B3). Our formulation suggests that terms which include second order derivatives should be 
considered on an equal footing with all other terms including a total of six derivatives, since 
they depend on the same number of fundamental quantum Green functions G^^ and . 

In a recent INM simulation of Ma and Stratt^^ the lowest order anharmonic contribution 
corresponding to Eq. (B6) was found to give a significant contribution to the total response 
of hquid Xenon. The fifth-order Raman response computed with this model did, however, 
still deviate significantly from that calculated using time correlation functions. ^^'^^ 

Both the diagonalization of the Hessian needed to obtain the normal modes and the eval- 
uation of the derivatives of the polarizability are time consuming processes. Which of these 
is the slowest depend on the model used to describe the polarizabihty. For the first order 
dipole induced dipole model the derivatives can be evaluated very effectively analytically. 
This simple model is, however, known to fail for systems with large polarizability^^ and in 
atomic liquids it gives an isotropic polarizability that is coordinate independent. 

In liquid phase, the basic assumption in the INM theory that motion can be described in 
static normal modes is not generally justified. In a liquid the molecules rotate and diffuse 
around changing the normal modes on the timescale that is studied in femto and picosecond 
experiments. Instead of using a static set of coordinates one can use a dynamic basis set. 
Ma and Stratt used a basis of normal modes that was allowed to change with time.^^ They 
obtained the Green's function describing the time evolution in one mode by applying WKB 



10 



theory. ^^'^^ In our notation, this result reads 

Mi.yuJi{T)uJi{0) 

The time evolution depends on the time dependent frequency u;j(r) or each normal mode. In 
this picture not only the frequencies are time dependent, but the polarizability derivatives 
and anharmonicities also change as the basis set change. Computationally it is expensive 
to repeatedly diagonalize the Hessian at short intervals to obtain the dynamic normal mode 
basis set and the method was not yet tested. 
Mode-coupling in k space: 

Reichman ^2, 53,54 Cao^^'^^ adopted Mode-Coupling theory of space nonlinear hydro- 
dynamics in k. Translational invariance then greatly simplifies the final expressions and the 
time dependence is accounted for through the dynamical structure factor F(k, t). They em- 
ployed the atomic first order dipole induced dipole model for the polarizabilit. The resulting 
classical expression is^^'^^ 



k 

X 



V^(k)V hT drs2 

1 rJTPfh- 1 

(28) 



kbT dT2i hT dT2i 

This result neglects the leading (fourth order) term (Eq. (B2)) and only retains one of 
the sixth order terms (Eq. (B3)). The time derivative of the dynamical structure factor is 
proportional to the Fourier transform of G'^~, while the dynamical structure factor itself 
is proportional to the Fourier transform of The V (k) / S (k)"^ factors correspond to 

polarizability derivatives in k space. ^^'^^ 

The quantum mechanical response function contains information about the local response 
assuming that the relevant coherence size is much smaller than the wavelength of the light. 
This is obviously the case for the response of many polyatomic molecules and aggregates. 
However, it usually holds for molecular liquids as well, where the coherence size underlying 
the response is small due to local disorder. It fails near critical points where the coherence 
lengths become very large. ^-^'^^''''^''''^ Since the response is local, this expansion should be 
made in real-space and not in k space. The simulation should provide a localized response 
since all six dipoles contributing to R^'^^^ should act within a small region (of order of the 
first solvation shell in CS2)) in order to generate the local response. Of course, the signal 
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should eventually be calculated in k space. This can be done purely macroscopically by 
solving Maxwell's equations within the local field approximation as described in chapter 16 
of reference 39. The only microscopic information that enters the signal is the effective local 
response, which can be obtained from simulations performed on the entire liquid. 
Mode-Coupling in real-space; Langevin equations: 

Using the Hamiltonian (Eqs. (7)- (10)) the Brownian oscillator motion is described by 
the following generalized Langevin equation.^^ 

MjQj{t) + MPjQ,{t) + MjY^ [ dT[^j,{t-r) + iEj,{t-T)U{T) = f j{t) + F,{t) (29) 

■ J-oo 

'-jij (Sjj) is the imaginary (real) parts of a self energy operator representing relaxation (level 
shift), fj is a Gaussian stochastic random force representing the bath degrees of freedom on 
the coordinate j and Fj is an external driving force. 

* a a 

E is related to 7 by the Kramers-Kronig relation. 

T.ij{u) ^ --Re r dJ^^^j^, (31) 



C"{uj) = Im (- 



In Appendix C the matrix of spectral densities is derived by solving Eq.(29).®°'^^ 

M{Q.^ + ujT.{uj) - lu^ + iuj^{uj))) ■ ^^^^ 
M, Vt and / are all diagonal matrices with matrix elements are given as follows Mij — SijMj, 
^ij — SijD,j and lij — 5ij. C"{uj) is the odd part of the spectral density, which is related 
to the even part C'{ijj) by the fluctuation-dissipation theorem. The matrix of spectral 
densities is given by C{ijj) = [1 -|- coth(/3^/2)]C"(a;), where 

^"^^^ ^ 7(Q^ - u^I + Euh-^iQ- - u-I + Eo;) + 7^0.^^^^"' ^^^^ 
Ordinary Langevin equations are obtained by first neglecting the frequency dependence 
of 7ij. Setting 7ij(a;) = 7ij and E.y(a;) = 0. Wc then take the overdamped limit ■yij >> flij 
of Eq. (32) where the matrix of spectral densities assumes the form^^ 

C"H = ^;;^AAc., (34) 

with the X matrices A and Q are A = 7-^^^ and A = M'^iQ-^f, ^^'^^ i^gj^g 
the number of Brownian oscillator modes. The nonlinear response of systems described by 
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Langevin equations can be obtained by using Eq. (34) for the matrix of spectral densities. 
Upon the substitution of Eq. (34) in Eq. (19) we get for the superoperator Greens functions 
in matrix notation 

G+-{r) = 2^(t) exp(-AT)AA (35) 

^ n=l 

where = 27m/hp are the Matsubara frequencies.'''^ In the high temperature {PhA « 1) 
hmit the second sum in Eq. (36) decays rapidly and may be neglected and the coth(/3/iA) 
factor can be approximated by 1/ f3hA in the first term, yielding 

G'++(t) = /3-iexp(-AT)AAA-^ (37) 

In the Generalized Langevin Equation approach of Kim and Keyes^^ the response function 
is factorized using the scheme suggested by van Zon and Schofield^^ which assumes fast 
decay of the fluctuating forces. This factorization results in the first order term for the 
2D response proportional to G~^~ {ts2)G^~ {t2i) (the second term in Eq.(B2)), whereas the 
first term G'"'"~(t32)G'+~(t3i) is neglected. The systematic mode-couphng factorization of the 
present paper is unambiguous and require no further assumptions about the behavior of the 
correlation functions. 

Adiabatic simulations: 

In systems with large scale variation of structure such as proteins and liquids it can be 
useful to employ a dynamic basis set instead of a static one. The natural starting point for 
such a description is to employ the adiabatic basis. ^^'^^'^^ A quantum description is obtained 
when one uses the eigenstates of the vibrational Hamiltonian as a dynamic basis. Using 
this basis we can derive an adiabatic theory for the response functions. When this basis 
changes slowly the adiabatic approximation can be employed, considerably simplifying the 
calculation of the time evolution, since the system remains in the same adiabatic state at 
all times. 

There are two independent Liouville space pathways contributing to the 2D response.^*^ 

R^^''\t32: T2l) = Ri{t^2, T21) + i?2(r32, T21) + C.C. (38) 

In the adiabatic approximation these are given by 

-Rl(T32,r2i) = (t) A^ca (t3 ) ha (^32 ) /^cb (^2) /fcq (t21 ) l^bajn) P ja) , 

^"■^ abc 
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R2{r32,r2i) = (^J ^IJ'cb{r3)hc{r32)lJ'cair2)Ibair2i)iJ,bairi)P{a), (39) 

^ ^ abc 

where 





rT2 




rT2 


Iab{r2i) = exp 


- / iea{T)dT 


exp 


/ ieb{T)dT 











eair) and ebir) are the eigenvalues of the Hamiltonian at time r. ^labi^) are the transition 
dipoles in the adiabatic basis. P{a) is the equihbrium population of state a. When the ener- 
gies change rapidly, the adiabatic approximation breaks down and nonadiabatic transitions 
need to be taken into account. 

In order to simulate the time evolution of the system one would have to diagonalize the 
vibrational Hamiltonian at short intervals. Such a treatment is therefore only feasible when 
the number of states is sufficiently small to allow repeated diagonalizations. This may be 
the case if a subset of vibrational coordinates are of interest and the rest are treated as bath 
coordinates. It should be noted that treating the vibrational states of the system explicitly 
only makes sense for high frequency modes, where the excited vibrational states have a low 
thermal population and only few states need to be considered. 

For all methods based on an expansion of the response function in the coordinates, 
evaluating the contributing terms gets increasingly more expensive with the order of the 
expansion. While the first and the second order derivatives of the polarizabihty needed to 
evaluate the lowest order terms can be handled rather easily, higher order derivatives get 
increasingly more difficult to evaluate. Unless the coordinates can be chosen such that the 
expansion may be truncated at some low order, it will be harder to simulate the classical 
response functions using the expansion, compared with real-space simulation methods. 



IV. CONCLUSIONS 



We have developed a systematic perturbative method for computing response functions 
using an expansion in the nonlincaritics, the effective dipoles, and the anharmonicities. 
Closed form expressions for the lower order terms have been derived for the first-, second- 
and third order response functions and applications were made to the ID and 2D response 
corresponding to third- and fifth order Raman techniques, respectively. 

The mode-coupling factorizations presented in this paper provide a unified framework for 
deriving all the approximate methods used so far in the simulations of the fifth order Raman. 
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Our expressions reduce to the instantaneous normal mode as well as the mode-coupling 
expressions by making additional approximations. This unified description allows a direct 
comparison of the various methods and can be used to develop semiclassical expansions. 

All existing simulations were compared and connected to the mode-coupling factorization 
presented in Sec. II. 

The present green function formalism allows the explicit incorporation of anhamonicities 
and gives a rigorous algorithm for truncating the expansion of the response functions, de- 
pending on the number of derivatives involved. At the same time, quantum corrections to the 
response functions may be computed as well. Mode-couphng theories in k space^^'^^'^^'^^'^^ 
utilize a non-local polarizability in order to describe response functions that are local in 
nature. This leads to the neglect of the leading order term in both the third and fifth order 
Raman response functions; contributions of all terms involving the first order derivative of 
the polarizability are missed. 

Collective coordinates are hkely to have Gaussian statistics. Identifying a set of collective 
coordinates that should allow a relatively simple interpretation of the response is the the 
key open challenge in the simulation of nonlinear response in the condensed phase. 
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APPENDIX A: SUPEROPERATOR ALGEBRA 

The following relations for the superoperators Q- and (5+ that follow directly from the 
definitions of the superoperators (Eq. (1)) can be used in expand the nonlinear response in 
terms of these superoperators. 



(Al) 
(A2) 
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(g^g^g% = ^(4Q^gigv + Q+Q-Q- + Q-W- + Q-Q-Q+) (A3) 

(g^g^gO- = {Q+QiQl + Q+QiQX + QtQ'+QX + ^QlQiQl) (A4) 

and so forth. Using these rules we can express the superoperator corresponding to an 
arbitrary product of ordinary operators as a product of superoperators 

where the coefficients /"*" and /~ are determined by apphcation of the rules (Eqs. A1-A3) 
and the i/'s denote either + or — . 

For a given set of expansion coefficients the n dimensional response can be expressed in 
terms of derivatives of the generating functional.^" 



(m + 2)!ni!---n„+i! V^. 

x,/"i) . . . T/™i ...V, 



tn+1 



X 



1 

9 



X 



^1 . . .^'"T- 



d d \ I d d 



9JkieS<) dJki^.emSO) \dJkTiTi.^'J ^Jki^^izS^L)^ 
■■■f^^...^m S{J{t)] (A6) 

APPENDIX B: 2D AND 3£) RESPONSE TO SIXTH ORDER IN THE PRIMARY 
COORDINATES 

The 2D and ?>D response functions are derived in the same way as the ID response 
function that was given in Eq. (24). The 2D response function expanded to sixth order in 
the primary coordinates is given by 

= i?gf ^ + i^g? + <^ + ^ + + • • • , (Bi) 

16 



where the first term is fourth order in the coordinates 

41? = E /^If /^i' VS'^G+- in,) (GJ- (rai) + Gj- (r^i)) . (B2) 

ijkl 

The remaining contributions are sixth order in the coordinates 

^2? = E 4Vi?/^^iG+-(r32)(G;-(r3i)G,i+(r2i) + G,V (B3) 



= ^ E f^f^&^n^G±^0){Gr-{T,,)Gt-{^^^^ 
ijklmn 

+ E (Gr(^32)G,V(r3i)G+^-(r32) +G++(r3i)G+-(r3i)G+-(r32)) 

+ E /^M/^i'UG'r(^32)Gr(r2i)G+-(r32)+Gr(r2i)G+-(^^^ 

ijklmn 

(B4) 

R?n = ^ E /^!>^Vi'^Gr(0)(G'r(^3i)G+-(r32) + Gr(r2i)G;:,^(^^^ (B5) 
4:nl = - E /^S'VJ'Vi'^^Si ^^I'^r (r3iOGi(ri'i)G+-(rr2) (B6) 

ijklmn ^ 

The two lowest-order terms contain numerous contributions. As noted earher^ aU terms 
in the 2D response contain nonhnearities either in form of anharmonicities or higher order 
derivatives of the dipole operator. 

The 3D response function expanded to sixth order in the coordinates is 

ijklmn 

+ G+-(r42)G;„-(r43)G+-(r2i) + G+-(T32)G+-(r2i)G+7(T« 

+ E /^S>S'V^Vi^^G',t-(r43)G;-(r42)G+-(r4i). (B7) 

ijklmn 

Numerous eight order terms exists including an anharmonic term. They can be readily 
obtained using the rules outlined earlier and will not be given here. 
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APPENDIX C: THE MATRIX OF SPECTRAL DENSITIES 



The matrix of spectral densities for a system described by the Hamiltonian defined in 
Eqs. (7), (10) and (11) can be determined by solving the generalized Langevin equation 
(Eq. (29).^^'^'^'^^ The generalized Langevin equation can also be written on the form 

M,g,(t) + M,Q,g,(t) + J2 dr^^I^^^^^^Q^ir) = ^ + F,{t) (CI) 

With the external driving force Fj{t) and the rapidly fluctuating force of the bath on the 
primary coordinate j fj{t). This force is determined from the system-bath Hamiltonian 
(Eq.(lO)). 

where the second term is independent of the bath coordinates. 

Comparing Eqs. (29) and (CI) allows us to identify 7 and E as the real and imaginary 
part of the correlations function of the fluctuating forces fj (t) . 

^ (f , V a .\ - {^fAt-r)Sfm) 

-V + -V = MkeT ^ 

Averaging over the bath coordinates and taking the Fourier transform of the generalized 
Langevin equation (Eq.(Cl)) gives: 

i 

In matrix form this gives: 

M(n2 _ ^2 J ^Y.uj- i-fu) {Q{uj)) = F{uj) (C5) 

The matrices M, and / are all diagonal with matrix elements Mjj = SijMj, ilij = Sijilj 
and lij = 6ij. Q and F are vectors. 

The change in the coordinates induced by an external driving force is 

(QM) = aMF(.) = ^^^._J^^_,^^f (-). (C6) 

which deflne the susceptibility a (a;). 

The odd part of the spectral density matrix is the imaginary part of the susceptibility: 

1 



C"{u;) = Im 
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(C7) 



where the imaginary part of a matrix is: Im A — (A — A^)/2i. 

The even part of the spectral density is related to the odd part by the fluctuation- 
dissipation theorem. 7 and E are determined by the correlation function of the fluctuating 
forces. The fluctuating force was determined by the derivative of the Hamiltonian. 

{Sfj{t)SMO)) = ^ cj^c,p{qa{t)qp{0)) (C8) 

The part of the fluctuating force that does not depend on the bath coordinates vanishes, 
when averaged over the bath coordinates. 

7 is determined by the real part of the time correlation function of the bath coordinates: 

^^^•(*) = E ^f^MQa{t)qpm (C9) 

a/3 ^ 



1 C jQ Cjp 

'^ijijjj) is the Fourier transform of 7ij(t) and is given in Eq. (30) 



y cos(a;„t) (CIO) 
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